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ABSTRACT 

We look at two contrasting spin-down models for isolated radio pulsars and, accounting 
for selection effects, synthesize observable populations. While our goal is to reproduce 
all of the observable characteristics, in this paper we pay particular attention to the 
form of the spin period vs. period derivative (P—P) diagram and its dependence on 
various pulsar properties. We analyse the initial spin period, the braking index, the 
magnetic field, various beaming models, as well as the pulsar's luminosity. In addition 
to considering the standard magnetic dipole model for pulsar spin-down, we also con- 
sider the recent hybrid model proposed by Contopoulos & Spitkovsky. The magnetic 
dipole model, however, does a better job of reproducing the observed pulsar popu- 
lation. We conclude that random alignment angles and period dependent luminosity 
distributions are essential to reproduce the observed P-P diagram. We also consider 
the time decay of alignment angles, and attempt to reconcile various models currently 
being studied. We conclude that, in order to account for recent evidence for the align- 
ment found by Weltevrede & Johnston, the braking torque on a neutron star should 
not depend strongly on the inclination. Our simulation code is publically available and 
includes a web-based interface to examine the results and make predictions for yields 
of current and future surveys. 
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1 INTRODUCTION 

The statistical properties of the radio pulsar population 
provide valuable constraints on the birth properties, 
evolution and radio lifetimes of neutron stars and 
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key diagram to reproduce is the logarithmic plot of pulsar 
periods, P, and their rates of slowdown, P, usually referred 
to as the P-P diagram. The many factors that affect the 
distribution of pulsars on the P—P diagram include the 
evolution of pulsars' spin-down, luminosity, age, magnetic 
field, beaming, and the angle between the magnetic field 
and the rotation axis. While each parameter individually 
affects how pulsars evolve across the diagram, covariances 
between multiple parameters also influence the results. 

The goal of this paper is to compare two pulsar spin- 
down models using Monte Carlo simulations of the Galactic 
young pulsar population. The first model we consider treats 
the pulsar as a rotating magnetic dipole from which there 



is a simple relat ionship between the brak ing torque and the 
magnetic field l|Goldlll968l : IPacinil Il968l ) . This model has 
been used ext ensively in pulsar population studies i n the 
past (see, e.g. [Lvne et al.lll985l : iNaravan &: Ostrikedll990l : 
iBhattacharva et al.lll992T) and has most r ecent ly been im- 
plemented by IFaucher-Giguere fc Kaspil (|2006t ). hereafter 
FK06. In this paper, we desire to undergo a direct com- 
parison of their results with our simulations. 

The other model we consider, by 

IContopoulos fc Spitkovskvl (|2006l b hereafter CS06, treats 
the pulsar as a combination of a misaligned magnetic 
dipole rotating in vacuum and an aligned magnetic dipole 
rotating in an atmosphere with ideal magnetohydrodynamic 
conditions. In addition to implementing these two models 
with parameters as described in detail in FK06 and CS06, 
we investigate the effects of changing one or more of the 
assumptions made by these authors in an attempt to obtain 
a more accurate interpretation of the relationships between 
these parameters and the resulting P-P diagram. 

As described in Section [2l we have created a com- 
puter simulation package which implements the above mod- 
els and allows us to simulate the distribution of pulsars in the 
Galaxy and change the individual properties of each one. By 
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"observing" this population with accurate models of com- 
pleted pulsar surveys (Section [3j, we are able to create a 
simulated sample and compare its properties to those of the 
real observed sample. In Section [4] we briefly describe our 
method for comparing two population samples. In Section 
we present and discuss our results. Our main conclusions 
are summarized in Section [6] To facilitate comparison with 
other studies, and provide a tool for the community, a ma- 
jor feature of this work is to make our simulation code freely 
available on a website which is discussed in Appendix A. 



2 MODELING THE UNDERLYING PULSAR 
POPULATION 

The computer program we have developed for this work, 
evolve, is an extension of the software package, psrpop 
which was initially developed to study the Galactic distri- 
bution of pulsars det ected in the Parkes Multibeam surveys 
|Lorimer et alj |2006). For that work, time-independent pul- 
sar populations were synthesized and compared to the ob- 
served sample. The main aim of this paper is to develop 
a time-dependent model of the population and is realized 
through a program we call evolve. In this section, we de- 
scribe the logical flow of the program which determines first 
whether a model pulsar is potentially observable (i.e. is ra- 
dio loud and beaming towards us), before computing any 
kinematic or spatial parameters. 

The first step is to randomly assign an age, t, and a mag- 
netic field, B, to each pulsar. The age is selected between 
and £ max , with t max being 10 9 years. As per FK06, the 
magnetic field is determined by a log-normal distribution, 
centred around a mean value, /ii og s = 12.65, with a stan- 
dard deviation a\ OE B = 0.55. The initial period, Po, is then 
randomly generated. Following FK06, Po is chosen from a 
Gaussian distribution centred around fip = 300 ms with a 
standard deviation ap = 150 ms. 

We next compute the inclination angle between the 
magnetic field axis and spin axis, X- I n implementations 
of the magnetic dipole model, it is usually assumed (see, 
e.g., FK06) that x — 90°. However, a more realistic model 
would account for a distribution of x an d use this informa- 
tion in a self-consistent way when calculating the beaming 
fraction discussed below. We explore a number of options, 
first fixing \ = 90° , and subsequently considering a random 
distribution of angles by selecting cos x from a flat distribu- 
tion between and 1. The latter choice effectively orients 
the magnetic axis as a uniform random vector over 4-7T sr. 
Finally, following recent em pirical evidence in favour of a 
secular decay of % with time (We ltevrede Sz Johnstonll2008l . 
hereafter WJ08), we also consider an exponentially decaying 
model defined as 



sinx = smxoexp(-t/t d ), 



(1) 



where xo is the initial angle (which we again choose by se- 
lecting cos xo from a 0-1 flat distribution), and td is the 1/e 
decay timescale. Note that we choose this decay law over 
the one considered by WJ08 for mathematical and compu- 
tational convenience. As shown in Fig. [T] the two depen- 
dencies give similar results. The braking index, n, is then 
assigned. In general, n is related to u, the frequency of ro- 
tation, v = 1/P, and its time derivative 




e(years) 



Figure 1. The time dependence of the alignment angle shown 
for Equation [TJ (dotted curve) and the model used in WJ08 (solid 
curve). We use 60° for a reference as roughly half of all angles will 
be above and half will be below that value for a given random 
distribution of angles. 



-Kv , 



(2) 



where K is a constant that de pends on the pulsar's mo ment 
of inertia and radius (see, e.g. lLorimer fc Kra mer 2005]). We 
can use a single index for all pulsars, or randomly assign an 
index to each pulsar. Specific choices depend on the spin- 
down model under consideration which are described below. 

Each spin-down model yields a time-dependent equa- 
tion for the pulsar's period and period derivative. The pul- 
sar's period evolves in time depending on the spin-down 
model, the magnetic field, the braking index, and the align- 
ment angle. Considering first the magnetic dipole model, the 
period dependence can be found by solving 



P n - 2 P = fcB 2 sin 2 X (£), 

where sinx(i) is defined by Equation [T] The constant 
8tt 2 P 6 



k = 



3/c 3 



(3) 



(4) 



where R is the radius of the neutron star, assumed to be 
10 km, and I is the moment of inertia, assumed to be 
10 38 kg m 2 . For this choice of parameters, k takes the nu- 
merical value 9.768 x 10" 40 , B is in units of Gauss while P, 
Po and t are in seconds. 

The implementation of this model in FK06 assumes n = 
3 and x = 90° . In this case, the period dependence 



P(t) = ^/pJT^kBH. 



(5) 



More generally, for any braking index other than n = 1, and 
a constant value of x> we find 



P(t) = [P n_1 + (n - l)kB 2 tsmx}~ 



(6) 



Finally, extending this result for a time dependent inclina- 
tion angle of the form assumed in Equation [T] we find 
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P(t) = 



+ 



(n- 



—t d kB 2 sin 2 xo ' 



-2t/t d \ 



■(7) 



Depending on the model under consideration, one of 
Equations [5H7] is used to calculate the pulsar's current pe- 
riod, P. Equation [3] is then used to find the correspond- 
ing value of P. At this point, a simple test is made to de- 
termine if the pulsar has crossed the so-called death line. 
The death line signifies when a pulsar becomes radio-quiet, 
and iBhattacharva et al.l (I1992T ) quantified it as the locus of 
points for which 

^- = 0.17 x 10 12 G s -2 . (8) 

If the pulsar has crossed the death line, the simulation marks 
the pulsar as dead (i.e. radio quiet) and moves on to cre- 
ate the next pulsar. If the pulsar has not crossed the line, 
then we continue on with the simulation. It should be noted, 
however, that a sm all number of pulsars (e.g. J2144— 3933; 
I Young et alii 19991 ) have crossed this theoretical line and are 
still observable at radio frequencies. 

For our implementation of the CS06 model, the tech- 
nique is slightly different. The time evolved period and 
period-derivative are calculated by integrating and then nu- 
merically solving 



P = 3.3 x 10" 



1 - 



p_\ 2 ' n / b \ 2 (-Py 1 
pj Vio 12 gJ \1s) 



■ cos 2 x 



where 



0.81 x 



(l0 12 g) (p ) 



(9) 



(10) 



In this case, the pulsar is declared dead when the period 
becomes greater than Pdeath- 

To account for the fact the pulsar radiation is only 
beamed to some fraction of 4ir sr, a variety of beaming mod- 
els have been imp lemented in evolve. The def ault model is 
from the work of iTauris fe Manchester! (1 19981 ) , henceforth 
TM98, who found an empirical relationship for the fraction 
of the whole sky illuminated by a pulsar: 



/(P) = 0.09 



loj 



(n) 



+ 0.03. 



(11) 



Other models from iBiggsl <fl990h. iLyne fc Manchester! 
l|l988l) , and iNaravan fc Vivekanand1(|l983h . henceforth B90. 
LM88, and NV83 respectively, determine the beaming frac- 
tion from the angular beam radius 



P = PoP 1 



(12) 



where the constants po and 7 are uniquely defined in each 
model. Assuming that radio waves are emitted from both 
poles of the pulsar and that the beams are circular beams, 
and following TM98 leads to a beaming fraction dependence 
on magnetic inclination angle of 



f(p,x) 



2 sin x sin p 
cos (x - p) 
1 - cos (x + p) 
1 



X > P;X + P < 2 

X>P,X + P>f 

x<p,x + p<| 

X< P,X + P> f ■ 



(13) 



model that assumes / = 0.2 for all values of P. To de- 
termine whether a pulsar is beaming towards us or not, we 
compare the value of / computed in one of the above ways 
with a random number between and 1. Only those pulsars 
for which this random number is less than or equal to / are 
deemed to be beaming towards us and we move on to the 
calculations of the final properties of the pulsar. 

To describe the location of model pulsars in the Galaxy, 
we use a regular Cartesian (x, y, z) coordinate system with 
the Galactic center at the origin and the position of the 
Sun at (0, 8.5, 0) kpc, again identical to the procedure per- 
formed by FK06. Each pulsar's initial x and y positions on 
the Galactic plane are determined by randomly choosing a 
point in the plane of our Galaxy along the spiral arms. We 
follow the procedure described by FK06 to implement the 
spiral arm structur e. For the radial d i stribu tion, we use the 
form described by lYusifov fc Kiicukl (2004). Then, disper- 
sion away from the plane is used to determine the initial 2 
position. For this, we use an exponential distribution with a 
mean 2o=50 pc, and then arbitrarily choose a sign. 

To model pulsar birth velocities, we use the optimal 
model found in FK06, in which the initial velocity com- 
ponents of each pulsar in the x, y and z directions are 
drawn from exponential distributions with a mean absolute 
value of 180 km s _1 . Although other distributions are avail- 
able within evolve, we do not investigate them any fur- 
ther here. To model the resulting evolution in the Galactic 
gravitational potenti al, we follow the method described by 
lLorimer et ail (|l993h . where the integration in the gravi- 
tat ional galactic potent i al is d one using a model described 
by ICarlberg fc Innanenl (|l987h . From the final position in 
the model galaxy, we compute the pulsar's distance from 
the Sun, d, and its expected dispersion measure (DM) and 
scatter-broadening timescale at 1 GHz. The latter two quan- 
tities are derived using the N E2001 model for the Galactic 
distribution of free electrons l|Cordes fc Lazioll2002f ). 

The radio luminosity is calculated, following FK06, with 
a P and P power law dependence, defined by 



log L = log L + a log P + 13 log (P/HT 15 ) + 8 L , 



(14) 



where Lq is 0.18 mjy kpc 2 , and 5l is randomly chosen from 
a normal distribution with as L =0.8. The program allows 
for a and /3, two free parameters, to be varied to determine 
the best functional dependence. 

For comparison, we also make available in the code a 
simple model in which the luminosity is independent of all 
other parameters. Following FK06, we take the probability, 
p(L), of a given luminosity as being 







L G [0 mjy kpc 2 , 0.1 mjy kpc 2 



Also available is a simple period-independent beaming 



p(L) oc ^ L"T5 L 6 [0.1 mjy kpc 2 , 2.0 mjy kpc 2 ) (15) 
L~ 2 L e [2.0 mjy kpc 2 , 00). 



3 MODELING THE OBSERVED PULSAR 
POPULATION 

The steps described above allow us to create a population 
of synthetic pulsars that are potentially observable. That is, 
they are deemed to be beaming towards the Earth and the 
spin-down models we are investigating predict that they are 
radio-loud. We wish to compare this population with those 
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pulsars that are potentially detectable by current surveys. 
Since the primary objective of this work is to reproduce and 
extend upon the work of FK06, in this paper, we focus on the 
sample of pulsars d etectable by the Parkes Multibeam (see 
lLorimer et alj|2006t and references t herein, hereafter PMB) 
and Swinburne ( Edwards et alJlioOll . hereafter SMB) multi- 
beam pulsar surveys. By accurately modeling the detection 
thresholds of these surveys as described below, and selecting 
only those model pulsars that are theoretically detectable, 
we can form samples of model observed pulsars. When FK06 
were carrying out their work, the sample of real pulsars de- 
tected in these surveys was 1065. In our analysis, we will 
use the sample of 1135 pulsars currently catalogued. Our 
criteria is to select all isolated pulsars not associated with 
a globular cluster and having both P < 10~ 12 and P > 30 
ms. 

To form our model observed samples, the first step is to 
compute the apparent flux d ensity of each mode l pulsa r, S. 
Following standard practice lLorimer fc Kramer! (|2005ri . we 
neglect geometrical factors in the inverse square law rela- 
tionship and find S from the pulsar's distance, d, and radio 
luminosity, L, as follows 

s =i- ^ 

Note that this flux density is denned (following our defini- 
tion of L) to be at an observing frequency of 1.4 GHz, the 
frequency at which the PMB and SMB surveys were carried 
out. For the purposes of this work, no assumptions about 
the pulsar spectral index distribution are necessary. 

To model the detection threshold of the pulsar surveys, 
it is also necessary to model the pulse widths. Following 
FK06, we assume that the intrinsic pulse width, Wi n t, is 
5% of the pulse period. The observed pulse width, W b B , 
and signal-to-noise ra tio, S/N, are th e n calc ulated using the 
method described in lLorimer et all l|2006h in their Equa- 
tions 1-7. Pulse w idth models from recent papers (see, e.g., 
ISmits et all 120091 ') have also been considered and are avail- 
able for use in the code. 

All pulsars which lie inside the sky boundaries covered 
by the two surveys, have S/N> 9, and where W b s < P 
are deemed detectable and saved for subsequent analysis, as 
described below. Our simulation proceeds until 1135 pulsars 
are detected, to match the combined sample found in the 
PMB and SMB surveys. 



4 METHOD FOR COMPARISON 

To compare our trial results quantitatively, we use the 
Kolmogorov- Smirnoff (KS) test, a non-parametric estima- 
tor based on the maximum deviation seen in the cumulative 
distr ibution of two sample data sets (see, e.g., iPress et al.l 
1986). For a given distribution, the KS test returns a statistic 
that is used to determine a probability, Q, that two samples 
came from the same underlying population. As described be- 
low, we make use of this probability to investigate whether 
model pulsar populations are inconsistent with the observed 
data. 

To investigate the sensitivity of the KS test to chang- 
ing model parameters, we run a single simulation using the 
optimal parameters from the FK06 model. We then gener- 
ate four further Monte Carlo realizations of the same model 



Table 1. Here we show how the KS probability for the period 
distribution behaves when given simulated samples that are ex- 
actly the same as, similar to, and completely different than an 
original sample. See text for details of the simulations. 



TRIAL 


Qp 


Qp 


Qp 


# 


(SAME) 


(SIMILAR) 


(DIFFERENT) 


1 


0.1299 


0.0296 


< 10~ 12 


2 


0.8928 


0.0094 


< io- 12 


3 


0.7397 


0.1098 


< 10" 12 


•1 


0.6958 


0.1744 


< IO" 12 


AVG 


0.6145 


0.0808 


< 10~ 12 



and compare them to the original simulation using the KS 
probability for the period distribution, Qp. Next, we make 
a small modification, changing \ip from 300 ms to 280 ms, 
run four more simulations and compute Qp by comparing 
with our original simulation. Finally, we make a more sig- 
nificant change to the model, changing to 250 ms, op 
to 100 ms, and HiogB to 12.00. We run four further simu- 
lations using these parameters and record Qp against the 
original simulation. From the results of these simulations 
that are summarized in Table [T] we conclude that individ- 
ual KS probabilities above 10% suggest that two popula- 
tions are consistent with another, while probabilities in the 
range 1-10% imply possible differences between populations. 
However, the steep decline in the KS probability when two 
populations differ significantly suggest that it is a powerful 
metric to weigh one model against another. 

In this work, we wish to track the multiple observed 
parameters relevant to our simulations. Therefore, following 
iBhattacharva et al.l (|l992l ). we define a figure of merit for 
each simulation 



FOM = log[Q P xQpxQixQ h 



(17) 



Here the subscripts on each KS Q value refer to the distribu- 
tions of period (P) , period derivative (P) , Galactic longitude 
(I) and Galactic latitude (b). 



5 RESULTS AND DISCUSSION 

In this section, we will describe and discuss our results. Our 
main goal is to compare the power law spin-down model 
(Eq. [2]) with the new spin-down model proposed by CS06 
(Eq. [9} under a variety of assumptions. The logical flow of 
our approach is summarized in Fig. 



5.1 Simulation Models 

We begin by using the basic model parameters from FK06 
and applying them to both spin-down laws which we call 
models 1A and IB. In this convention "A" refers to the 
power-law spin-down model and "B" refers to the CS06 
spin-down model. Model 2 attempts to improve the origi- 
nal model by changing the luminosity law and using random 
inclination angles to better match the observed sample. Sub- 
sequent models are motivated to be more physically realistic 
with the incorporation of: angle- dependent beaming models 
and a range of braking indices for models 3 A and 3B, and 
inclination angle decay for models 4 A and 4B. 
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Figure 2. Schematic showing the logical flow of our approach to comparing the power law (A) and CS06 spin-down (B) models. We begin 
in model 1 by using the best parameters from FK06. As described in the text, subsequent models then contain additional complexity in 
an attempt to be more physically plausible representations of the true population. 



Table [2] provides an overall summary of our results, 
showing base- 10 logarithms of the individual KS probabili- 
ties and the FOM for each model. We also compute the mean 
birth rate (72.) of pulsars required to produce the observed 
sample sizes for each model. This is defined as the ratio of 
the total number of pulsars generated in each simulation to 
the maximum age of the population, t max . 

5.1.1 Model 1: FK06 basic parameters 

Adopting the optimal parameters found by FK06, we gen- 
erate model A which is shown in Fig. [3] The histograms 
show a relatively accurate replication of the original FK06 
results. While we are generally able to reproduce the results 
of FK06, we find somewhat lower KS probabilities for P and 
P than FK06. Obtaining an exact match to results obtained 
with a different simulation is challenging and the discrepan- 
cies we find highlight the difference in the implementation of 
the model in the two codes. We duplicate most of the major 
results well. Some other parameters, such as a and /3, have 
to be changed slightly from the values presented in FK06, 
but we obtain the majority of their results using our code. 
These changes are addressed in the remaining models. 

For model B, we retain all the parameters of model A 
with the exception of the spin-down law, which we replace 
by our implementation of CS06 described in Section [2] As 
evidenced by the results shown in Table [2] the P and P 
distributions are not as well reproduced. 

This model has mixed results. Many areas, such as the 
effect of the alignment angle, the effect of having no angle 
dependence, and the underlying distribution of pulsars lying 
below the death line (see Fig. 1, 3, 6 of their paper), can be 




Galactic Longitude [DegJ Galactic Latitude [Deg] 



Figure 3. Histograms of some of the key properties from model 
1A. The solid histograms represent the observed pulsar popula- 
tion. The dotted histograms are the results of our simulation. 
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Figure 4. Our resulting P—P diagrams, with the adopted death 
line shown in each plot for reference. Shown in the upper panels 
are our implementations of FK06 (model 1A, upper left) and CS06 
(model IB, upper right). The lower panels show our improved 
simulation (model 2A, lower left), which is based on an FK06 
spin-down model, and the real observed sample of pulsars (lower 
right). 

replicated with our code. However, the KS statistics are not 
as high as they are for the FK06 spin-down model. This 
applies not only to Model 1, but to the rest of the models 
as well. 



5.1.2 Model 2: Modified luminosity law and random 
inclination angles 

In an attempt to improve the low KS statistics seen in the 
P and P distributions for model 1, we search a range of 
different luminosity indices. Using Equation 1141 we vary a 
from —1.9 to —0.8 and /3 from 0.1 to 0.9. Our best result 
occurs when a and /3 take the values of —1.0 and 0.5, while 
FK06 found optimal values of —1.5 and 0.5. As discussed in 
FK06, simulation techniques vary from model to model, and 
small deviations may occur between models in the search for 
optimal parameters. 

So far we have assumed orthogonal rotators in the spin- 
down models. A more realistic and self-consistent approach 
is now followed by assigning random inclination angles to all 
pulsars. 

The results from model 2 are significantly better than 
model 1, as seen in Table [2] At this point, we can see how 
much improvement is gained by changing only two param- 
eters, so we move on to test other parameters to see if this 
improvement continues. 

5.1.3 Model 3: Beaming model and braking index 
variations 

Motivated by the improvements seen in model 2, we now 
relax the requirement for the braking index n = 3 and allow 
a range of braking indices. Based on the observed sample, a 
sensible choice is to select n from a fiat distribution in the 
range 2.5 < n < 3.0. 

Our previous models adopt the TM98 beaming model, 
which is used in FK06. However, we now consider computing 



the beaming fraction based on a period-dependent opening 
angle law. 

These results are worse than the previous model in 
terms of KS statistics, however, they are still better than 
the original model. For model A, the birthrate decreases 
from nearly 4 pulsars per century to a more acceptable 3 
per century. However, model B has just the opposite trend: 
the birthrate changes from 3 to 4 pulsars per century. 



5.1.4 Model 4-' Inclination angle decay 

Our final model is motivated by the recent results of WJ08, 
who found evidence for alignment of the spin and magnetic 
axes. In both models 4A and 4B, the initially random incli- 
nation angles decay according to Eq. [T]with an exponential 
decay timescale of td = 10 r yr. Model 4A uses the spin- 
down law found in Eq. [Jj while model 4B assumes the CS06 
relation given in Eq. [9] As shown in Table 2, both models 
perform very poorly, with significantly lower figures of merit 
than the earlier models without alignment. For the case of 
model 4A, this poor agreement can be understood when it 
is seen that Eq.[jjis essentially equivalent to an exponential 
decay of the magnetic field — more generally, an exponen- 
tial decay of the braking torque. From the work of FK06, 
and preliminary simulations we have carried out, we know 
that a decaying braking torque is inconsistent with the ob- 
servations. For the case of model 4B, the strong dependence 
on the spin-down model (Eq. [9j) with \ is also the reason for 
the low figure of merit. 

These results present a dilemma. On one hand, the em- 
pirical evidence for alignment presented by WJ08 appears 
to be very robust. On the other hand, we have found that 
magnetic dipole spin-down laws provide a good description 
of the population without inclination angle decay. The spin- 
down laws we have implemented depend critically on this 
inclination angle. This implies that one way to resolve the 
discrepancy is to modify the spin-down law so that inclina- 
tion angles are removed. While there is no physical basis for 
taking the sin 2 \ term out of Eq. [3] we find that a modified 
version of model 4A which uses Eq. [5] instead of Eq. [Jj does 
produced improved results, though the overall FOM (-6) is 
still lower than seen for the other models. In a future paper, 
we intend to investigate this issue further. Based on our cur- 
rent results, however, it appears that the braking torque on 
a pulsar is independent of its magnetic inclination angle. 

Finally, for completeness, we have compared models 
2A, 3A, and 4A w ith the pulse width model mentioned in 
ISmits et al.l (|2009T ). The birthrates for all three trials remain 
consistent with each model, however, the KS statistic FOM 
is generally worse. We conclude that changing the way in 
which pulse widths are modeled will not improve our re- 
sults. 



5.2 Individual parameters 

Separate from our optimal models, we test various pulsar 
parameters by themselves in order to gain a better under- 
standing of how they affect the resulting pulsar population. 
In this section, we take a closer look at each of these param- 
eters. 
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Table 2. Summary of KS statistics, figures of merit (FOM) and birth rate (TV) comparing FK06 with the four models considered in 
this paper. See text for details of the various models. 



Model 1 Model 2 Model 3 Model 4 



Parameter 


FK06 


A 


B 


A 


B 


A 


B 


A 


B 


P 


-2.15 


-1.68 


-12.0 


-0.94 


-7.00 


-1.75 


-12.0 


-12.0 


-10.7 


P 


-1.20 


-4.41 


-5.35 


-0.51 


-0.67 


-0.08 


-1.76 


-6.37 


-2.21 


I 


-0.90 


-0.12 


-0.12 


-0.07 


-0.09 


-0.28 


-0.40 


-0.10 


-0.25 


b 


-1.14 


-0.50 


-0.13 


-0.23 


-0.30 


-0.95 


-1.87 


-1.01 


-1.56 


FOM 


-5.39 


-6.71 


-17.6 


-1.75 


-8.06 


-3.06 


-16.0 


-19.5 


-14.7 


Ti (psrs/century) 


2.8 


3.85 


1.97 


3.56 


2.84 


2.97 


3.75 


5.31 


3.68 



5.2.1 Beaming 

We run multiple simulations to test the beaming models. 
For each trial, we alter only which beaming model is to be 
tested, and leave all other parameters the same. The con- 
stant beaming fraction model, or one in which 20% of all 
pulsars beam towards Earth, and the NV83 beaming model 
both give poor results. The KS statistics using these models 
are very close to zero for all simulations. 

The beaming models of LM88, B90, and TM98 each give 
much better results. The TM98 model yields the best pop- 
ulation of pulsars, but it is not statistically better than the 
other two models. Finally, we try a period-dependent model, 
similar to one described in WJ08. This gives us our best re- 
sult overall result in terms of KS statistics and birthrates. 

Based on these results, we can say a requirement for the 
population is a period-dependent beaming model, as well 
as small opening angles, similar to LM88, B90, and TM98. 
These three beaming models, along with WJ08, are all sta- 
tistically about the same. 



5.2.2 Braking Index 

As a test of braking indices, we use our simulation Model 
2A. We try indices ranging from 3.2 to 2.5, as well as a very 
realistic random model in which indices are chosen at ran- 
dom between 2.5 and 3.0. We find very few differences when 
varying the braking indices, as shown in Table [3] Due to 
the lack of any significant difference between various brak- 
ing indices, we can adopt the realistic random model for use 
in our simulations. 

It is useful to again look back at the A models from 
Table [2] We notice a decrease in birthrate when going from 
model 2 to model 3. The two parameters that are altered be- 
tween the two models are the braking index and the beaming 
model. The results presented in Table |3] make it clear that 
this drop in birthrate is not entirely due to the braking in- 
dices, and thus the angle-dependent beaming model must 
affect the birth rate as well. 



5.2.3 Luminosity 

We run many luminosity simulations by changing the power 
law indices of P and P. The best result we obtain yields 
indices of —1.0 for P and 0.5 for P. Our trials also include 
a simulation involving random luminosities, i.e. no period 
dependence. The results are less than ideal, and not close to 



Table 3. We show our results from 9 simulations with different 
values of the pulsars' braking indices. This is an implementation 
of the our Model 2A, changing only the braking index during each 
simulation. The resulting birth rate, in units of pulsars born per 
century, and FOM (as described previously) is shown for compar- 
ison purposes. 



Braking Index Ti (psrs/century) FOM 



2.5 < n 



3.2 


3.75 


-3.78 


3.1 


3.99 


-4.09 


3.0 


3.56 


-1.75 


2.9 


3.51 


-1.63 


2.8 


3.52 


-2.27 


2.7 


3.53 


-3.66 


2.6 


3.25 


-1.41 


2.5 


3.28 


-2.74 


< 3.0 


3.60 


-1.38 



the results we obtain when using a period-dependent lumi- 
nosity distribution. Thus, we conclude that there must be 
some period dependency for the luminosity. Similar conclu- 
sions were reached by FK06, however, their optimal value 
for the P power law index was —1.5. 

For our ideal simulations, our underlying luminosities 
have a distribution that is very similar to the log-normal 
distribution found in Fig. 15 of FK06. In this regard, we 
note that this distribution provides, in our opinion, the best 
current estimate of the pulsar lumin osity function. Unlik e 
previous power-law models (see, e.g.. lLorimer et all [2006), 
the log-normal does note require a lower luminosity cutoff 
and is strongly recommended for studies which require some 
reasonable assumption about the luminosity function. 

5.2.4 Covariances 

With such a large number of model parameters in our simu- 
lations, we are mindful during our study to keep track of any 
interdependencies that might be expected. The biggest co- 
variances we find involve the luminosity. Changing the brak- 
ing index along with the luminosity function can be done in 
such a way that two very similar results can be obtained with 
two completely different luminosity distributions and brak- 
ing indices. For example, n=3.0 and power indices of —1.5 
and 0.5 could yield the same results as n=2.7 and indices 
of —1.2 and 0.8. A similar relationship is found between the 
luminosity and the alignment angle. Running a simulation 
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with all angles equal to 90° might give similar results to 
another simulation with all angles equal to 45° and having 
luminosity indices of —1.8 and 0.9. Most of these covariances 
are beyond the scope of this paper, simply due to the large 
matrix of parameters that would be analysed. 

5.2.5 Death line 

For completeness, we run simulations of both models that 
do not use a death line. As one might expect, we obtain 
an overabundance of high period pulsars that have very low 
luminosities. Our results here just confirm that, within the 
framework of the models we have developed here, a death 
line is indeed required when modeling a distribution of pul- 
sars. While it might be possible to construct a model in 
future which does not require a death line, we are not cur- 
rently aware of any straightforward means of achieving this. 



6 CONCLUSIONS 

We have successfully implemented two pulsar spin-down 
models and accounted for selection effects as far as possi- 
ble to synthesize populations of observable radio pulsars. 
From a statistical comparisons of various models, we have 
gained a better understanding of how their parameters af- 
fect the spin-down evolution of pulsars on a P—P diagram. 
Our main conclusions are as follows: 

• The magnetic dipole spin-down model from FK06 works 
best with our simulations. We are able to reproduce the 
results from their paper and obtain some improved results 
by modifying some of their parameters. The model of CS06 
produces poorer results. Further modifications of this model 
appear to be required to improve it. 

• The braking index, n, does not have a significant impact 
on our results. Having a braking index of n=3.0 works just 
as well as the rather unphysical scenario of n >3.0. For our 
optimal simulation, we use the most realistic model, which 
picks a random braking index between 2.5 and 3.0 for each 
pulsar. 

• The optimal configuration for magnetic inclination an- 
gles is a random distribution. Models in which the incli- 
nation angle decays on timescales of ~ 10 7 yr do not re- 
produce the observations well. Further modifications to the 
spin-down laws appear to be required in order to account for 
the strong empirical evidence for alignment found by WJ08. 

• Pulsar luminosities must have a period dependence. We 
are able to use a power law for the period and period deriva- 
tive to replicate the P—P diagram. By tweaking the expo- 
nents in that law, we can alter the resulting distribution of 
pulsars on our diagram. Even with period dependent lumi- 
nosity laws and beaming models, a death line is required to 
explain the dearth of pulsars in the lower right of the P-P 
diagram. 

Studying these relationships between the P-P diagram and 
the various pulsar parameters allowed us to become more 
aware of covariances, eliminate a few parameter models, and 
overall obtain better insight to the behaviour of pulsars as 
they evolve across the P—P diagram. Future studies will en- 
able us to further narrow down some parameter models and, 



in particular, allow us to investigate the issue of magnetic 
alignment. 
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APPENDIX A: AN OPEN-SOURCE 
APPROACH TO PULSAR POPULATION 
SYNTHESIS 

Follow ing the initial version described by lLorimer et al.l 
(2006), the source code used in this work is freely avail- 
able at http://psrpop.sourceforge.net We briefly de- 
scribe the new features of the pulsar simulation package, 
psrpop, as well as announce the availability of a website, 
http : //psrpop . phys . wvu . edu , used to investigate past and 
future pulsar surveys. 

The biggest new feature of the software package is the 
ability to evolve a pulsar's spin period in time. Currently, 
spin-down models from FK06 and CS06 are available for use 
in the simulations. By having this evolution feature, we can 
now watch a pulsar's life cycle, from birth to death, and ob- 
serve how its period changes with time. The program also 
has the ability to use various beaming models, luminosity 
laws, and alignment angle functions. These allow for further 
studying of the individual pulsars, and a more complete un- 
derstanding of pulsar populations as a whole. 

The psrpop website currently has the capability of sur- 
veying any of the eight model populations generated for this 
paper, and future models will also be available for use. The 
user can "search" these model populations using previous 
surveys such as the Parkes Multibeam Survey, theoretical 
surveys such as one using the proposed Square Kilometer 
Array, or they can create their own survey. One of the main 
benefits of running a custom survey is the ability to predict 
yields of future surveys. 

After surveying one of the pulsar populations, the web- 
site outputs the total number of pulsars detected, a P-P 
diagram, and and few comparison plots. These plots contain 
histograms of the pulsar properties detected in the survey 
that overlay histograms of the observed pulsar population. 
Some possible properties that can be shown in the histogram 
plots are pulse period, period derivative, dispersion measure, 
and flux density. 



